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ABSTRACT 

We present an implementation of smoothed particle hydrodynamics (SPH) with im¬ 
proved accuracy for simulations of galaxies and the large-scale structure. In particular, 
we implement and test a vast majority of SPH improvement in the developer version 
of GADGET-3. We use the Wendland kernel functions, a particle wake-up time-step 
limiting mechanism and a time-dependent scheme for artificial viscosity including high- 
order gradient computation and shear flow limiter. Additionally, we include a novel 
prescription for time-dependent artificial conduction, which corrects for gravitation¬ 
ally induced pressure gradients and improves the SPH performance in capturing the 
development of gas-dynamical instabilities. 

We extensively test our new implementation in a wide range of hydrodynami- 
cal standard tests including weak and strong shocks as well as shear flows, turbu¬ 
lent spectra, gas mixing, hydrostatic equilibria and self-gravitating gas clouds. We 
jointly employ all modifications; however, when necessary we study the performance 
of individual code modules. We approximate hydrodynamical states more accurately 
and with significantly less noise than standard GADGET-SPH. Furthermore, the new 
implementation promotes the mixing of entropy between different fluid phases, also 
within cosmological simulations. 

Finally, we study the performance of the hydrodynamical solver in the context of 
radiative galaxy formation and non-radiative galaxy cluster formation. We find galac¬ 
tic disks to be colder and more extended and galaxy clusters showing entropy cores 
instead of steadily declining entropy profiles. In summary, we demonstrate that our 
improved SPH implementation overcomes most of the undesirable limitations of stan¬ 
dard GADGET-SPH, thus becoming the core of an efficient code for large cosmological 
simulations. 

Key words: hydrodynamics - methods: numerical 


1 INTRODUCTION 

Smoothed particle hydrodynamics (SPH) is a commonly em¬ 
ployed numerical method in astrophysics. It solves the fluid 
equations (Landau & Lifshitz 1959) in a Lagrangian mass- 
discretised fashion, which ensures Galilean invariance and 
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conservation of mass, momentum, angular momentum, en¬ 
ergy and entropy. It was pioneered by Gingold & Monaghan 
(1977) and Lucy (1977) and has since then become one 
of the cornerstones of computational astrophysics. The dis¬ 
cretisation of mass automatically adapts spatial resolution 
by removing the constraint of handling geometry explic¬ 
itly. It also easily couples to N-Body schemes for calcula¬ 
tion of gravitational forces (Hernquist & Katz 1989). An 
excessive amount of papers and literature about SPH has 
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been produced over the past decades. We point out the lat¬ 
est reviews by Rosswog (2009), Springel (2010b), Monaghan 
(2012) and Price (2012a) for the basic concepts and e.g. 
Ritchie & Thomas (2001) for an extension to multi-phase 
fluids and Rosswog (2015) for a special-relativistic adaption. 
As every numerical method, SPH comes with its own set of 
benefits and pitfalls, which we address in this paper. 

The inability of traditional SPH methods to treat con¬ 
tact discontinuities and to mix different fluid phases is a 
long-standing problem (Agertz et al. 2007; Wadsley et al. 
2008). It leads to a completely numerical spurious sur¬ 
face tension at the discontinuities preventing particle move¬ 
ments. Consequently, it results in a failure of these for¬ 
mulations of SPH to resolve fluid instabilities such as the 
Kelvin-Helmholtz or Rayleigh-Taylor instabilities (see e.g. 
Junk et al. 2010; Valcke et al. 2010; McNally et al. 2012; 
Puri & Ramachandran 2014). In applications to cosmic 
structure formation it causes entropy profiles to diverge to¬ 
wards the centres of dark matter haloes, at variance with 
Eulerian codes that predict entropy plateaus to build up 
(Frenk et al. 1999; Wadsley et al. 2008; Planelles & Quilis 
2009; Vazza 2011; Power et al. 2014; Biffi & Valdarnini 
2015; Rasia et al. 2015). This difference is due to the lack 
of mixing in simple SPH, which makes low-entropy gas in 
merging substructures sink toward the centre of the main 
structure. Many modifications have been proposed to over¬ 
come this problem. For example, Wadsley et al. (2008) pro¬ 
pose a mixing solution, which resolves the differences in the 
entropy profiles of dark matter haloes between Eulerian and 
SPH codes (Frenk et al. 1999). Further cosmological appli¬ 
cations have been performed by Shen et al. (2010). Firstly, 
the equation of motion (EoM) can be re-formulated from 
a standard ’density’ approach into a ’pressure’ based ap¬ 
proach (Saitoh & Makino 2013; Hopkins 2013). While the 
’pressure’ formulation correctly treats contact discontinu¬ 
ities, it leads to increased noise at strong shocks. Sec¬ 
ondly, considerable effort has been made to unite grid-based 
solvers for the fluid equations with the Lagrangian nature 
of SPH. Eulerian Godunov methods (see e.g. Cha et al. 
2010; Springel 2010a; Murante et al. 2011) and their cou¬ 
pling to Lagrangian methods is a promising alternative. 
Connecting a Lagrangian moving-mesh with grid-based 
solvers (Springel 2010a) or mesh-free approaches (Hopkins 
2015; Hopkins & Raives 2015) represent more advanced ap¬ 
proaches. Thirdly, artificially modelled conduction (AC) of 
internal energy can be employed to overcome the mixing 
problem. Most modern SPH codes include AC of some sort 
to diffuse entropy (Read & Hayfield 2012) or energy (Price 
2008) across particles. The use of AC has to be taken care¬ 
fully and it is only desirable at contact discontinuities in 
traditional ’density’ SPH and at shocks in modern ’pres¬ 
sure’ SPH. The application of AC in other regions can have 
catastrophic impact on the fluid dynamics and can smear out 
gravitationally established pressure gradients, thus leading 
to totally numerically induced transport of internal energy 
(Valdarnini 2012). 

Next, traditional SPH has difficulties treating subsonic 
turbulence as it experiences a high effective viscosity, which 
limits the inertial range (see e.g. Bauer & Springel 2012). 
Thus, traditional SPH cannot achieve high Reynolds’ num¬ 
bers compared to, for exmaple, Eulerian methods. This high 
effective viscosity is a function of resolution, but no general 


solution has yet been proposed to resolve this issue in gen¬ 
eral. For the correct capturing of shocks numerically mo¬ 
tivated artificial viscosity (AV) is commonly employed. It 
smooths the particle velocity distribution and gives order 
to the fluid sampling. However, AV is only desired at the 
shock and the fluid should be inviscid otherwise. Too much 
AV smears out physical motions and damps subsonic and 
turbulent motions in isolated tests (Bauer & Springel 2012) 
as well as in cosmological simulations (Dolag et al. 2005). 
Therefore, several different implementations of AV reduction 
are proposed (Morris & Monaghan 1997; Cullen & Dehnen 
2010). They are all based around a proper shock detection 
method and a time-dependent viscosity decay scheme. Ap¬ 
plication of such advanced schemes give better results in the 
description of fluid motions (Dolag et al. 2005; Price 2012b). 

Finally, it might seem easy to simply reduce the quan¬ 
titative errors by increasing the number of neighbours, 
which contribute to the local density and force estima¬ 
tors. However, the standard weighting functions of SPH re¬ 
spond differently to an increase of smoothing neighbours 
and can possibly become unstable to the pairing instabil¬ 
ity (Schuessler & Schmitt 1981; Price 2012a). Therefore, re¬ 
cently, alternative kernel functions immune against this in¬ 
stability are proposed for better fluid sampling and conver¬ 
gence (Read et al. 2010; Dehnen & Aly 2012). The advan¬ 
tage of flexible geometry of SPH comes with difficulties in 
creating well-defined initial conditions or sampling analyti¬ 
cal profiles, where we use either glass set-ups (White 1996) 
or Weighted Voronoi-Tesselations (Diehl et al. 2012). 

To overcome the named disadvantages we implement a 
large set of improvements for SPH into the developer ver¬ 
sion of the cosmological N-Body / SPH simulation code 
GADGET-3 (Springel et al. 2001; Springel 2005). We in¬ 
clude a time-step limiter for strong shocks, a time-dependent 
viscosity scheme for subsonic turbulence, a high-order gra¬ 
dient estimator and shear flow limiter for shearing motions, 
an improved kernel function for convergence and a time- 
dependent artificial conduction scheme to promote fluid mix¬ 
ing. We discuss the accuracy and the performance of our new 
scheme in hydrodynamical standard test problems, within 
quiet and violent environments as well as in Idealized simu¬ 
lations of galaxy and galaxy cluster formation, in which our 
new scheme is applied to reasonable astrophysical problems. 

The paper is organised as follows. The improved imple¬ 
mentation of hydrodynamics is presented in Section 2. In 
Section 3 we validate our SPH algorithms in a set of hydro- 
dynamical standard tests and we proceed to standard tests 
with gravity in Section 4. We continue in Section 5 with 
Idealized applications to the evolution of an isolated disk 
galaxy and a forming galaxy cluster. We summarise our de¬ 
velopments and code performance in Section 6. 


2 A NEW SPH IMPLEMENTATION 

We start with a presentation of the main equations cor¬ 
responding to a ‘standard’ and our ‘new’ formalism of 
GADGET-SPH. The formalism of SPH is already well de¬ 
scribed by a large number of reviews (see e.g. Price 2012a). 
We refer to the ‘standard’ version of SPH as the implemen¬ 
tation within the GADGET-3 code without our modifica¬ 
tions. We point out our modifications and discuss the ker- 
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nel function, the EoM, the particle wake-up scheme and the 
time-dependent AV and AC. 


2.1 Original code platform 

We implement our SPH modifications into the developer 
version of the cosmological N-Body/SPH code GADGET-3 
(Springel et al. 2001; Springel 2005). We evolve entropy as 
the thermodynamical variable (Springel & Hernquist 2002) 
and use the prescriptions for radiative cooling, supernova 
feedback and star formation following Springel & Hernquist 
(2003). In the following sections we compare two different 
SPH schemes (see also Table 1), which are distinguished as 
follows. The ‘standard’ implementation corresponds to the 
developer version of GADGET-3 without our modifications 
(Springel 2005). The ‘new’ implementation includes all the 
SPH improvements presented in this section. In principle, we 
always use the entire new scheme and we also employ the 
same set of all numerical parameters throughout our entire 
simulation test suite. Thus, unless otherwise stated, we do 
not tune individual standard tests or astrophysical applica¬ 
tions. However, if necessary we sometimes switch off some 
of the modifications to analyse their individual and isolated 
impact on several of the test problems. 


2.2 Kernel functions and density estimate 

Foremost, there is the question in a Lagrangian method how 
to derive fluid field quantities from a given set of point 
masses. In particular, the estimation of the gas density is 
crucial as many further equations rely on it. We employ 
the standard estimator of SPH and calculate the density 
p(x;) = pi of an individual particle i at the position x; by 
summing the contributions of N% neighbouring particles j 
within a smoothing radius /i(x,) at a distance Xij in a mass- 
weighted (rrij) and distance-weighted (Wij(xij, hi)) fashion 

Pi = m j Wij(xij,hi). (1) 

3 

Simultaneously, the smoothing length hi is a function of den¬ 
sity 


/ \ i/3 

ft(xi) = V , (2) 

where p defines the ratio of smoothing length to the mean 
distance between particles. Eqs. (1) and (2) roughly en¬ 
sure constant mass resolution throughout the simulation 
and have to be solved in parallel. This mimics the evolu¬ 
tion of spheres of the same mass 4/37T hi/pi = Nirrii but with 
varying number of neighbours. The number of neighbours 
varies across space and time with an increase or decrease of 
smoothing length and local quality of fluid sampling by the 
point masses. The weighting function is commonly chosen to 
decrease monotonically with distance, yield smooth deriva¬ 
tives, is symmetric with respect to Xij = Xji and has a flat 
central portion. A historical choice (Monaghan & Lattanzio 
1985) of kernel function 

Wij(xij,hi) =w(q)/hi, (3) 


is the cubic B-spline function with q = Xij/hi and 

of l-6q 2 + 6q 3 0 <q<\ 

w(q) = - l 2 (1 - qf | < q < 1 , (4) 

^ { 0 1 <q 

which we commonly employ with a choice of 64 neighbours 
in three dimensions. However, this traditional kernel func¬ 
tion is subject to the pairing (or clumping) instability when 
the number of neighbours is too large (see Price 2012a). An 
alternative choice to achieve better numerical convergence is 
necessary and an entire new family of kernels is needed. In 
a kernel stability analysis Dchnen & Aly (2012) show that 
the Wendland kernel functions are a much better choice. We 
choose the Wendland C 4 (WC4) kernel with 200 neighbours 
in three dimensions as our smoothing function without the 
pairing instability problem. The functional form of the C 4 
is given by 

w(?) = ^(!-g) 6 ( 1 + 6 9+y9 2 )- (5) 

For values of q > 1 it is set to w(q) = 0. The Wendland 
functions require similar computational effort as the cubic 
spline kernel but nevertheless, the total computational time 
increases due to larger number of neighbours. Therefore, we 
do not employ the higher-order C 6 functions because of the 
required 295 neighbours. In summary, the total computa¬ 
tional cost of the density and hydrodynamical force calcula¬ 
tion increases by a factor of about two in comparison to a 
cubic spline with 64 neighbours. However, a better estimate 
of the kernel will result in a more accurate density estimate 
and improved gradient estimators. These estimators are the 
cornerstones of the SPH formalism and determine the accu¬ 
racy and convergence rate in all our test simulations. 


2.3 Equation of motion 

The EoM for a system of point masses are derived (see 
e.g. Price 2012a) from a discretised version of the fluid La¬ 
grangian 



( 6 ) 


where v denotes velocities and u internal energy of individ¬ 
ual particles. The Lagrangian nature of SPH, when com¬ 
plemented with a symplectic time integration scheme, au¬ 
tomatically conserves mass, momentum, angular momen¬ 
tum, energy and entropy. We use the standard kick-drift-kick 
Leapfrog time integration of GADGET (Springel 2005). The 
EoM then follows from the principle of least action, where 
the spatial derivative of internal energy comes (if constant 
entropy is assumed) from the first law of thermodynamics 
dU = —PdV. We choose a volume element depending on 
density (V = m/p) and an adiabatic equation of state for 
the pressure P = Ap 1 , which is defined individually for ev¬ 
ery particle. We integrate entropy A as the thermodynamical 
variable of choice and thus employ what is commonly called 
’density-entropy’ SPH. 

The EoM in the ’density-entropy’ (for a derivation see 
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Standard 

New 


Density estimator 
Kernel function 
Neighbours (3D) 

Traditional 

Cubic spline 

64 

Bias-corrected 
Wendland C 4 

200 

Dehnen & Aly (2012) 

Dehnen & Aly (2012) 

Dehnen & Aly (2012) 

Equation of motion 
Grad -h terms 

Density-Entropy 

Yes 

Density-Entropy 

Yes 

Springel & Hernquist (2002) 

Springel & Hernquist (2002) 

Velocity gradients 
Artificial viscosity 
Balsara limiter 

Low-order 

Standard (constant) 
Low-order 

High-order 
Adaptive (locally) 
High-order 

Price (2012a); Hu et al. (2014) 

Dolag et al. (2005); Cullen &; Dehnen (2010) 

Balsara (1995); Cullen & Dehnen (2010); Price (2012a) 

Artificial conduction 
Hydrostatic correction 

No 

No 

Adaptive (locally) 
Adaptive (locally) 

Wadsley et al. (2008); Price (2008) 

Price (2008); Valdarnini (2012) 

Particle wake-up 

No 

yes (fu, = 3) 

Saitoh & Makino (2009); Pakmor (2010); Pakmor et al. (2012) 


Table 1 . Comparison of the ‘standard’ (column 2) and ‘new’ (column 3) SPH implementations in the GADGET code. Furthermore, we 
give some references (column 4) for extended descriptions and discussions. 


Springel & Hernquist 2002) for the hydrodynamical force 
part of an individual particle reads 


= ~y^j 

fi°*ViWij(hi) + fT^iWij{hj) 

hy<l j 

V Pi Pj \ 


(7) 


where the factor / co is a correction factor, which accounts 
for the mutual co-dependence of smoothing length h(p) and 
density p(h) and their corresponding derivatives. Its func¬ 
tional form is given by 


hi dpi 

3pi dhi 


( 8 ) 


Eq. (7) leads to a non-vanishing force at contact disconti¬ 
nuities even when pressure is constant. This is the artificial 
’surface tension’ of SPH, which suppresses particle move¬ 
ment across contact discontinuities. In the following sections, 
we present our equations in notation of internal energy u, 
which is related to the entropic function A = (7 — l)u/p 1 ^ 1 . 


2.4 Artificial viscosity 

2-4-1 Smoothing of jumps 

By construction, SPH solves the ideal Euler equation and 
no dissipative terms are included but those are necessary 
to describe discontinuities correctly. In highly dynamical re¬ 
gions (e.g. in shocks) fast particles commonly penetrate into 
regions of resting particles causing unwanted particle disor¬ 
der and oscillations in the sampling of the fluid. However, 
SPH already contains an intrinsic remeshing force but to 
re-establish particle order and capture shocks properly an 
additional dissipative term in velocity is needed. This AV 
aims to remove post-shock oscillations and noise and helps to 
smooth the velocity field (see Monaghan & Gingold 1983). 
We include AV in an energy conserving way with a contri¬ 
bution to the EoM of the form 

-tt = g ~ “ Vi ) v a Fi i’ ( 9 ) 

vise j 


and with a contribution to the energy equation of the form 


dm 

dt 


1 \ - m-j . ,2 v 


f shear sig, Vi 
Jij u ij 1 


( 10 ) 


where the symmetrised variables represent pij — (pi + pj)/2 
for the density, a^j = (cr” + a ])/2 as a numerical coef¬ 
ficient to include AV (see below) and /fj 168,1 = (/j ?hear + 
/| hear )/2 as the Balsara (1995) shear flow limiter (see Sec¬ 
tion 2.3.2 below), which aims to ensure the application of 
AV only in strong shocks (high velocity divergence) and not 
in rotating or shearing flows (high velocity curl). Further¬ 
more, in the above equation Fij = ( Fij(hi ) + Fij(hj))/2 
is the symmetrised scalar part of the kernel gradient terms 
ViWij(h-i) = FijVij/nj , which are used to linearly inter¬ 
polate the second-order Laplacian derivative in the velocity 
field diffusion equation. The pairwise signal velocity u^ B,v 
(first introduced by Monaghan 1997, and already used in 
GADGET-2) determines the strength of AV and directly in¬ 
cludes a quantitative measure of particle disorder 


l#^=C?+C$-ftx«, (11) 

where c s is the sound speed of the particles and p-ij = v;, • 
x y /xij with a commonly chosen pre-factor of ft = 3. AV 
is only applied between approaching pairs of particles (i.e. 
p-ij < 0) and otherwise switched off. The local signal velocity 
v“ s (also used by the time-step criterion, see Section 2.6) 
represents the maximum value of u“- s,v between all particle 
pairs ij within the entire smoothing sphere of particle i. 

The calculation of the viscosity coefficient af is based 
on an approach developed by Cullen & Dehnen (2010) but 
modified for more efficient computation as follows. The pres¬ 
ence of a shock is indicated via computation of velocity di¬ 
vergence contributions across the entire smoothing length 
by 

Hi — — Y, sign(V ■ v)jmjWij, (12) 

pi j 

where a shock corresponds to Ri ~ —1. In principle, an accu¬ 
rate calculation of Ri for every particle requires the previous 
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computation of (V-v)i for every particle. Therefore, an extra 
SPH summation loop added between the calculation of den¬ 
sity (where velocity divergence can also be calculated) and 
hydro forces would be necessary. For computational reasons 
we use the velocity divergence calculated in the previous 
time-step. Furthermore, a convergent flow is also indicated 
by a high velocity divergence but that condition does not dis¬ 
tinguish between pre-shock and post-shock regions. There¬ 
fore, we employ the time derivative of velocity divergence to 
determine a directional shock indicator 


Ai = £imax(0,-(V ■ v)i), (13) 

which is able to distinguish between pre-shock and post¬ 
shock regions. We calculate (V ■ v)i via interpolation be¬ 
tween the current and the previous time-step (as suggested 
by Cullen & Dehnen 2010) in the time interval A ti. 

Subsequently, we use the shock indicator Ri to deter¬ 
mine the ratio of strength of shock and strength of shear 
in quadratic form via 


|2(l-flQ 4 (V-v)i| 2 n4 x 

|2(1 — f?i) 4 (V ■ v)i| 2 + |V x v|? ’ 1 ’ 

which is proposed by Cullen & Dehnen (2010) as an addi¬ 
tional limiting factor for AV in Eq. (13) and was experi¬ 
mentally determined. Now, for every particle we can define 
and set the target value a l ° c,v of AV with the help of the 
directional shock indicator to 


loc.u 

OLa — Qimax 


h lAi 


h fAi + « 


i si SY 


(15) 


In the case, where the viscosity coefficient a/ is smaller than 
a* oc, ’\ we set the coefficient to a° c ’ v . Otherwise, we let it 
decay with time according to 


/ loc.i; v \ 


(16) 


which we integrate in time together with the hydrodynami- 
cal quantities. The constant l specifies the decay length and 
in our test problems we find a numerical value of l = 4.0 to 
give reasonable results. 


2.4-. 2 Gradient estimators 

We use the Balsara (1995) form of the shear viscosity limiter 


/f 


|V-v|j 

|V ■ v|j + |V x v |i + at ’ 


(17) 


with ai = O.OOOlCi /hi for numerical stability reasons. At 
this point, the question arises how to calculate the diver¬ 
gence and vorticity. The common curl estimator of SPH 
reads 


(V x v)j = —— Ym, (vj — Vi) x ViWij, (18) 

pi i 

which takes the lowest order error term into account. Since 
higher-order error terms are neglected, this formation per¬ 
forms very poorly in the regime of strong shear flows. 
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Therefore, we resort to a higher-order calculation scheme 
of the velocity gradient matrix (see similar approaches by 
Cullen & Dehnen 2010; Price 2012a; Hu et al. 2014). We fol¬ 
low the approach presented in Price (2012a) for the compu¬ 
tation of the gradient matrix. We expand Vj for every vector 
component k in a Taylor-series around i with 

Vj = Vi + ( dsVi)(xj - Xi) s + 0(h 2 ). (19) 

Inserting Eq. (19) into Eq. (18) yields an easy solution for 
the linear term dsvlf and the velocity gradient matrix by 
solving the following system of equations: 


X*' 9 = ~ x i )“VfW ij 


( 20 ) 


dvf 

dx a 


E m j( v i - Vi) fc VfWij, 


( 21 ) 


which requires a matrix inversion for ■ Conveniently, the 
estimator is also independent of density and thus, can be 
calculated in the same computational loop along with densi¬ 
ties. Subsequently, the updated estimates of velocity diver¬ 
gence and curl are calculated directly from the full velocity 
gradient matrix via 


(V-v)i 


dvr 

dx a ’ 


( 22 ) 


(Vxv)‘ = e^A. (23) 

In our test problems we find the low-order estimator of ve¬ 
locity divergence to give already satisfying results (see also 
appendix A2 in Schaye et al. 2015). In contrast, the low- 
order estimator of velocity curl performs very poorly and 
we obtain significantly improved results with the high-order 
curl estimator of Eq. (23). The high-order estimators are 
not restricted to the AV scheme but they also enter various 
other modules of the code, where their precise calculation is 
required. For example, this additionally greatly improves the 
approximation of fluid vorticity written into the simulation 
snapshots. 


2.5 Artificial conductivity with gravity correction 

2.5.1 Smoothing of jumps 

We move on to address the mixing problem in SPH by 
introducing a kernel-scale exchange term for internal en¬ 
ergy transport. We include AC for purely numerical rea¬ 
sons to treat discontinuities in the internal energy (similar 
to the capturing of velocity jumps by AV), which arise from 
our ’density-entropy’ formulation of SPH. We note that a 
’pressure-entropy’ formulation of the EoM is also able to 
address the mixing problem but it also requires the pres¬ 
ence of AC in order to smooth noise in internal energy 
behind shocks (Hopkins 2013; Hu et al. 2014). Thus, in ei¬ 
ther flavor of SPH the inclusion of AC is recommended and 
many different formulations of the AC diffusion equation 
have been investigated so far. Although their precise details 
vary across the literature, they all ensure conservation of in¬ 
ternal energy within the kernel. Price (2008), Price (2012a) 
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and Valdarnini (2012) propose the diffusion of internal en¬ 
ergy, while Read & Hayfield (2012) propose the diffusion of 
entropy. Wadsley et al. (2008) propose a first mixing formu¬ 
lation to resolve the differences in entropy profiles within 
cosmological comparison simulations (Frenk et al. 1999) be¬ 
tween grid and SPH codes. The diffusion coefficient is ap¬ 
proximately proportional to a c v sie,c Xij and the numerical 
coefficient a c is commonly treated as constant through space 
and time. We adapt the formulation of a spatially varying 
coefficient of Tricco & Price (2013) and additionally calcu¬ 
late a limiter depending on the local liydrodynamical and 
gravitational states. We compute the gradient of internal 
energy as 


(Vu)i = — ^2 m A u i — Ui)ViWij (24) 

pi j 

and approximate the AC coefficient 


hi |Vu|i 

3 |u»| 


(25) 


as a measure of noise of internal energy sampling on kernel 
scale. The time evolution (i.e. spatially varying SPH dis¬ 
cretization of the second-order diffusion equation) of the in¬ 
ternal energy for each particle and its neighbours is then 
given by 


a; )n) j( :y- c (26) 

j Pij 

where we employ the choice of Price (2008) for signal velocity 
depending on the pressure gradient of the form 


dui 

dt 


sie.c 
17. . 


I Pi - A 


Pij 


(27) 


and ay = (a£ + aJ)/2 is the symmetrised conduction 
coefficient, which are individually limited to the interval 
[0,1], In the literature several other forms of AC (see eg. 
Wadsley et al. 2008; Valdarnini 2012) or approaches to the 
mixing problem (see e.g. Hopkins 2013) have been proposed. 


2.5.2 Gravity limiter 

We note that the amount of AC applied depends on the gra¬ 
dients of internal energy and of pressure. In the case that 
the thermal pressure gradient is determined by gravitational 
forces (i.e. hydrostatic equilibrium) this method would in¬ 
correctly lead to unwanted conduction. In the following, we 
determine the contribution of hydrostatic equilibrium to the 
total thermal pressure gradient and present a method to 
limit the amount of conduction. Firstly, for every individual 
active particle, we project the gravitational force F 9 onto 
the hydrodynamical force F^ 1 and calculate the partial force 
Ff of Fj 1 , which is balanced by Ff to 


F? = 


(Ff'F' 1 ) h 
|F?| 2 i ' 


(28) 


The sign of Ff depends on the spatial orientation of the 
force vectors. Secondly, we subtract/add the partial force 


F 9 from/to the hydrodynamical force F(“ and obtain F 9 , 
which we call the gravitationally adjusted hydrodynamical 
force 

F; = F^ 1 + F 9 (29) 

and which we use to determine a limitation factor 8f for AC 



(30) 


The limiter ensures that AC is only applied to the part of 
F^ 1 which is not balanced by F 9 . The exponent q represents 
a scaling for the aggressivity of the gravity correction. We 
limit our correction factor to the interval [0, 1] and directly 
multiply it onto the individual AC coefficients of. The lim¬ 
iter performs only as well as the hydrodynamical scheme is 
able to resolve hydrostatic equilibrium (in the ideal case the 
angle between force vectors is 180°). However, in SPH sim¬ 
ulations small-scale noise is present at all times within the 
kernel and thus also in the force vector angles. The expo¬ 
nent q (applied after the boundary verification) can then be 
understood to account for the noise in the particle distri¬ 
bution and mimics an opening angle of force vectors. After 
extensive studies and performing a variety of test problems, 
we settle with q = 5. The limiter returns zero in the case 
no hydrodynamical forces are present and one in the case, 
where no gravitational forces are present. We are aware that 
in the presence of strong pressure gradients and rotational 
forces our approach only marginally limits the amount of 
AC applied. However, we did not encounter major problems 
in our simulations performed with the ‘new’ scheme so far. 
Therefore, we assume this issue to be not too important at 
the present state. 


2.6 Particle wake-up and time-step limiter 

GADGET employs individual time-steps for all of the parti¬ 
cles to increase computational efficiency. Thus, the particle 
population is split into a set of active particles, whose hydro¬ 
dynamical properties are integrated in the current time-step 
and a set of inactive particles, which reside on larger time- 
steps. These individual time-steps are computed from the 
local thermodynamical properties of each particle. However, 
the splitting between active and inactive computational re¬ 
gions creates problems, where both sets of particles are over¬ 
lapping. In the case of a rapid gain in velocity or entropy an 
active particle can penetrate into a region of inactive parti¬ 
cles. The inactive particles do not notice the sudden presence 
of the highly dynamical particle and therefore large gradi¬ 
ents in the time-steps and unphysical results can occur. As 
a treatment we adopt a time-step limiting particle wake- 
up scheme as proposed and implemented in the GADGET- 
3 code by Pakmor (2010) and Pakmor et al. (2012) with 
the help of K. Dolag. It is similar to the time-step limiting 
scheme described by Durier & Dalla Vecchia (2012) and can 
be considered an extension of the Saitoh & Makino (2009) 
mechanism. Furthermore, our limiter compares signal veloc¬ 
ities instead of time-steps and accounts for incorrect extrap¬ 
olations. For every active particle, in every time-step, the 
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individual time-steps themselves are re-computed according 
to 


where C is the Courant factor and i>; lg the maximum signal 
velocity (see Section 2.4). For the calculation of the time- 
step the maximum of the signal velocity computed between 
the active particle i and all its neighbour particles j within 
the entire kernel is used. GADGET employs a check during 
the hydrodynamical force computation for large differences 
in signal velocities (see Eq. (11)) within the kernel by eval¬ 
uating 

vlf > Uvf e , (32) 

with a tolerance factor f w corresponding to a wake-up trig¬ 
gering criterion, which captures sudden changes in the pair¬ 
wise signal velocity. From our hydrodynamical standard 
tests we find f w = 3 to give reasonable results. Additionally, 
the fluid quantities of the recently woken-up particles could 
have already been predicted half a time-step into the future. 
Therefore, the incorrect extrapolation is removed and the 
contribution from the real time-step added. These correc¬ 
tions are performed for all particles for which the time-steps 
are adapted. 


3 HYDRODYNAMICAL TESTS WITHOUT 

GRAVITY 

We evaluate the performance and accuracy of the two differ¬ 
ent SPH implementations with a first set of standard prob¬ 
lems. These first test problems are purely hydrodynamical 
and do not include gravity or more advanced physics, yet. 
Throughout all the test problems we use an adiabatic in¬ 
dex of 7 = 5/3, the same set of numerical parameters (see 
Section 2) and we do not specifically tune individual test 
problems. 

3.1 Sod shock tube 

We consider the Sod shock tube problem (Sod 1978) to study 
the SPH behaviour in a simple weak shock test. We set up 
630000 particles of equal masses using a relaxed glass Hie in a 
three-dimensional periodic box with dimensions Ax = 140, 
Ay = 1 and A z = 1. On the left half side of the computa¬ 
tional domain (x < 70) we initialize 560000 particles with a 
density of pL = 1.0 and a pressure of Pl = 1.0. On the right 
half side of the computational domain (x > 70) we initialize 
70000 particles with a density of pR = 0.125 and a pressure 
of Pr = 0.1. 

Fig. 1 shows the results of the test problem at time 
t = 5.0. In general, both SPH schemes agree fairly well with 
a reference solution (green line) obtained with the ATHENA 
code (Stone et al. 2008) but we note the following differ¬ 
ences. In the ‘standard’ scheme (blue dots), the disconti¬ 
nuity in internal energy results in a ’blip’ of pressure (see 
Fig. 2) and energy, which corresponds to the artificial spu¬ 
rious surface tension of SPH. The issue of the ’blip’ has 
been discussed for a long time (see e.g. Monaghan 1992). In 


the ‘new’ scheme (red dots), AC promotes mixing, resolves 
the discontinuity, regularizes the pressure and provides a 
treatment of the ’blip’. A closer look at individual particles 
(see Fig. 3) shows that the pure noise in velocity of parti¬ 
cles behind the shock front is lower, which is a direct result 
of the improved prescription of AV. However, reducing the 
viscosity gives rise to post-shock ringing. Additionally, the 
change of kernel function improves the sampling quality of 
the fluid and yields a smoother estimation of density. At 
last, the time-step limiter is of little importance due to the 
weak shock in this test. However, as seen in Fig. 2 the ’blip’ 
is not completely removed and this is where some residual 
surface tension shows up. 


3.2 Sedov blast 

We consider the Sedov blast problem (Sedov 1959) to study 
the SPH behaviour in a simple ultrasonic strong shock test. 
We set up 130 3 particles of equal masses using a relaxed 
glass file in a three-dimensional periodic box with dimen¬ 
sions Ax = Ay = Az = 6 kpc. In the entire computational 
domain we initialize the particles with a density of 1.24 x 10 6 
Mq kpc -3 and one Kelvin as temperature. At the centre of 
the box we point-like distribute the energy E = 6.78 x 10 53 
erg to mimic a supernova explosion among the nearest 102 
particles. 

Fig. 4 shows thin slices through the centre of the simula¬ 
tion box and Fig. 5 the corresponding particle distribution 
at time t = 0.03. Furthermore, we perform the ‘standard’ 
scheme test also with the time-step limiter ( f w = 8000) be¬ 
cause of the very strong shock (Mach 3> 100) of the blast 
and otherwise any comparison will fail. Without the lim¬ 
iter, shocked particles penetrate into quiescent regions caus¬ 
ing a highly distorted fluid sampling, which results in an 
incorrect solution leading to an incorrect propagation of 
the shock front (see discussions in Saitoh & Makino 2009; 
Durier & Dalla Vecchia 2012) and corresponding smooth¬ 
ing. The entire ’new scheme’ reproduces the analytical so¬ 
lution (black line) very well, with the ‘standard’ run (green 
dots) totally failing, and the ‘new’ run (red dots) captur¬ 
ing the position, density and temperature of the shock fairly 
well. In addition, we show a partially improved ‘standard’ 
run (blue dots), where we enabled the time-step limiter but 
nothing else. This run also yields reasonable good results 
in this test, but as we see later comes short in other tests. 
We see that the ‘new’ scheme yields a smooth distribution 
of particles within the central region and therefore a well- 
resolved, but smoothed due to AC, temperature solution. 


3.3 Keplerian ring 

We consider the Keplerian ring problem (Cartwright et al. 
2009; Cullen & Dehnen 2010) to study the SPH behaviour in 
a simple rotating and shearing test problem. We set up 20000 
particles of equal masses sampling a two-dimensional ring 
with a Gaussian surface density profile with a peak at radius 
R = 15.0 kpc and a standard deviation of a = 2.0 kpc. For 
numerical reasons we initialize the distribution in concentric 
shifted circles and not in a random fashion (Cartwright et al. 
2009). We set the particles on Keplerian orbits around a 
central 10 9 Mq point mass with a rotation period of t = 2n. 
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Figure 1. Sod shock tube. We show the spatial distribution 
of particles (every 10 th particle is plotted) for density, thermal 
pressure, total energy and velocities at time t = 5.0. Both SPH 
schemes capture the shock well but with differences as follows. 
The ‘new’ scheme converges better in the density estimate and 
the presence of AC nearly removes the pressure blip at the contact 
discontinuity. 
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Figure 2. Sod shock tube. We show a zoom-in on the pressure 
blip (see Fig. 1, upper right panel) at time t = 5.0. Only a small 
residual of surface tension is left. 
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Figure 3. Sod shock tube. We show a zoom-in on velocity in re¬ 
direction (see Fig. 1, bottom right panel) at time t = 5.0. Besides 
some post-shock ringing, the post-shock noise is smaller with the 
’new’ scheme. 
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Figure 4. Sedov blast wave. We show thin slices through the 
centre of the computational volume at time t = 0.03 of the test 
performed with the ‘new’ scheme. The shock front is clearly visible 
in the gas density (left panel) as well as the temperature (right 
panel). 




Radius [kpc] 

Figure 5. Sedov blast wave. We show the radial distribution of 
particles (every 5 th particle is plotted) at time t = 0.03 with a 
time-step limiting criterion of f w = 3. We have performed the 
‘standard’ run with a time-step limiting criterion of f w = 8000 
(green lines, otherwise no meaningful comparison can be per¬ 
formed) and also f w = 3. The classic ‘standard’ scheme (green 
lines) fails to capture the shock, while the ‘new’ scheme (red lines) 
captures the position and evolution of the shock front much bet¬ 
ter compared to the analytical solution (black lines). This test 
shows the importance of the time-step limiter. 
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We choose the sound speed orders of magnitudes smaller 
than the orbital velocity to ensure thermal stability of the 
ring. In contrast to our common set of numerical parameters, 
we start without a minimal value of AV because it would 
immediately trigger instability. In the absence of AV the o 
ring should be stable. ‘g 

Fig. 6 shows the results of the test problem at the times 
of onset of runaway instability. We perform all test runs with 
the WC4 kernel in order to exclude possible effects caused 
by the smoothing scale of kernel sizes and differential esti¬ 
mators. Due to the highly sub-sonic nature and the absence 
of strong shocks, the impact of AC and the time-step limiter 
is negligible. The initially stable ring (top left panel) evolves 
as follows for different implementations of AV. In the ‘stan- c 
dard’ scheme (top right panel), the ring is only stable for 
about two dynamical times, before the instability has fully 
developed and the ring breaks up. Also, the Balsara limiter 
does not succeed in limiting AV because of the insufficient 
calculation of vorticity. In the ’M&M’ scheme (bottom left 
panel) we use the implementation of a low-viscosity scheme 
initially proposed by Morris & Monaghan (1997) and imple¬ 
mented into GADGET by Dolag et al. (2005). Their scheme 
uses a time-dependent evolution of numerical AV coefficient 
Q.'l to suppress AV in the absence of shocks and manages to 
keep the ring stable for about seven dynamical times. How¬ 
ever, the M&M scheme requires a minimum value of AV and 
also uses a low-order estimator for vorticity, which leads to 
the ring break-up. At last, we show the results of the ‘new’ 
scheme (bottom right panel), which we used without a min¬ 
imum value for AV. This scheme uses a high-order estimator 
of velocity gradient matrix, which results in a very accurate 
calculation of divergence and vorticity. Therefore, also the 
computation of the Balsara shear flow limiter is very accu¬ 
rate and suppresses AV completely within the entire ring 
structure. We do not note an artificially induced transport 
of angular momentum and orbital changes of test particles. 
Consequently, the ring remains stable for many dynamical 
times and the initial Gaussian surface density distribution 
is preserved until we stopped the simulation. 


3.4 Cold blob test 

We consider the blob test (Agertz et al. 2007; Read et al. 
2010 ) set up with publicly available initial conditions 1 to 
study the SPH behaviour in a test problem with interacting 
gas phases and surfaces. We initialize 9641651 particles of 
equal masses using a relaxed glass file in a three-dimensional 
periodic box with dimensions Ax = 10, Ay = 10 and Az = 
30 in units of the cold cloud radius. A cold cloud is centred 
at x, y, z = 5 and travels at a Mach number of M = 2.7. The 
background medium is set-up ten times less dense and ten 
times hotter than the cloud. Spherical harmonics are used to 
seed large-scale perturbations onto the surface of the cloud. 
Because of the low Mach number shocks we expect the time- 
step limiter to be only of minor importance. 

Fig. 7 shows thin slices through the density structure at 
various times. In the ‘standard’ scheme, the cold gas cloud is 
prevented from dissociating by the follow major effect (see 


1 http://www.astrosim.net/code/doku.php 
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Figure 6. Keplerian Ring. We show the initial set-up (top left 
panel) as well as the results of three different AV methods. For a 
fair comparison we only vary the AV scheme and none of the other 
SPH modifications. The importance of AV becomes clear as high 
amounts of viscosity lead to numerical accretion of particles onto 
the central point mass. Because angular momentum is conserved, 
the ring breaks up and an instability develops. With a low-order 
Balsara limiter, neither the standard SPH viscosity (top right 
panel) nor the M&M viscosity (bottom left panel) are able to 
preserve the ring. The ‘new’ scheme (bottom right panel), which 
has a the time-dependent AV coupled with a high-order limiter, 
is able to preserve the ring to even very late times. 


also e.g. Agertz et al. 2007). The presence of artificial sur¬ 
face tension confines the blob of cold gas. This is clearly 
visible by the numerically induced stretching of the cloud. 
Cold material, which should have been mixed into the am¬ 
bient hot medium is confined within an elongated structure. 
In the ‘new’ scheme, the presence of AC promotes the mix¬ 
ing with external ambient medium. We also note that the 
‘new’ scheme resolves different shock structures propagating 
through the box both more accurately and smoothly. 

Fig. 8 illustrates the dissipation of the gas cloud by 
tracking the time evolution of cold blob mass. We associate 
(see also Agertz et al. 2007) particles with the cold cloud 
with a temperature criterion of T < 0.9 • Text (in contrast 
to the external ambient medium) and a density criterion 
of p > 0.64 • pci (in contrast to the initial density of the 
cloud). In the ‘standard’ scheme (blue line), only half of the 
cold gas mass is mixed into the hot ambient medium over 
five dynamical times. The effects of the improved AV and 
WC4 kernel (green line) are negligible. The major impact 
and contribution to cloud dissociation is made by AC (pink 
line) and the corresponding introduced mixing process. The 
results are close to a test run performed with the ENZO 
(Bryan et al. 2014) code in a comparable set-up, which we 
took from Hopkins (2013). However, some residual surface 
tension remains. 
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Figure 7. Blob test. We show thin slices of gas density through the centre of the computational domain at times t = 0.0, 3.0, 6.0 and 
10.0. In the ‘standard’ scheme, numerical surface tension prevents mixing between cold and hot phases leading to an artificial stretching 
of the cloud and an unphysical solution. In the ‘new’ scheme, AC helps to promote cloud dissociation; however, some residual surface 
tension remains left. Furthermore, the shock structures throughout the box are more defined and better resolved. 


3.5 Kelvin-Helmholtz instability 

We consider the Kelvin-Helmholtz instability (Agertz et al. 
2007; Read et al. 2010) from publicly available initial con¬ 
ditions 1 to study the SPH behaviour in a simple shearing 
instability test. We set up 1548288 particles of equal masses 
using a cubic lattice in a three-dimensional periodic box with 
dimensions Ax = 256 = Ay = 256 and Az = 16 kpc, which 
is centred around (0,0,0). In the central half of the box 
(|y| < 64) we initialize 512000 particles with a density of 
pi = 6.26 • 10 3 MQkpc -3 , temperature of Ti = 2.5 • 10 6 K 
and a velocity in ^-direction of Vi = —40 km/s. In the outer 
half of the box (|y| > 64) we initialize 1036288 particles 
with a density of p 2 = 3.13 • 10 2 Mgkpc -3 , temperature of 
T 2 = 5.0 • 10 6 K and a velocity in ^-direction of V 2 = 40 
km/s. To trigger the instability, we perturb the velocity in 
y-direction with a mode of wavelength 128 kpc and ampli¬ 
tude of 4 km/s at the boundary layer that is exponentially 
damped towards the upper and lower edge of the box. 

Fig. 9 shows thin projections through the specific en¬ 
tropy structure of the test problem at various times. In the 
‘standard’ scheme, the fluid evolves in a laminar fashion and 
the growth of perturbations is totally suppressed by the ar¬ 
tificial surface tension confining the central gas stream and 
large amounts of AV damping velocity perturbations (see 
also e.g. Agertz et al. 2007; Price 2008). In the ‘new’ scheme, 


the high-order Balsara shear limiter successfully limits AV 
and allows large-scale perturbations to develop two promi¬ 
nent roll-ups. Additionally, AC nearly removes the artifi¬ 
cial surface tension between the two gas phases and pro¬ 
motes mixing within the roll-ups. The entire test set-up does 
not evolve completely symmetric because of small secondary 
perturbations caused by the initial set-up on a cubic lattice. 
Most importantly, the high-order AV and AC prove crucial 
for this test problem, while the WC4 kernel and time-step 
limiter are of less importance. At the late stages, in this 
set-up the ‘new’ scheme is dominated by diffusion. 

3.6 Decaying Subsonic Turbulence 

Recent comparisons of standard SPH implementations with 
static and moving mesh grid codes have sparked a debate 
about the capabilities of SPH to model subsonic turbulences 
(see e.g. Bauer & Springel 2012; Price 2012b; Hopkins 2013, 
2015). We study the behaviour of our ‘new’ scheme in Ide¬ 
alized simulations of decaying subsonic turbulence. In par¬ 
ticular, we are interested in the effective viscosity of the two 
schemes and the behaviour of the ’SPH noise’ under con¬ 
ditions appropriate to galaxy formation and cluster simula¬ 
tions, i.e. non-isothermal, decaying motions from solenoidal 
and compressive modes. As most baryons on cosmological 
scales are in weakly collisional plasmas, numerical models 
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Figure 8. Blob test. We show the fraction of cold gas as a func¬ 
tion of time for four different SPH schemes. In the ‘standard’ 
scheme, numerical surface tension prevents mixing between the 
cold and hot phases. The importance of AC becomes clear as it 
promotes mixing between gas phases, which allows a dissociation 
of the cloud. Only some residual surface tension is left. Compar¬ 
ing to a run performed with the ENZO grid code, we find our AC 
scheme to model mixing in a conservative way. 


should aim to minimize viscosity where possible (see e.g. 
Brunetti & Lazarian 2007). 

3.6.1 Grid and particle conversion procedures 

We set up 512 3 particles of equal masses within a periodic 
box of side length 3 Mpc/h using carefully relaxed SPH 
glass files to minimize spurious initial kinetic energy. Sub¬ 
sequently, we define a velocity field on a grid of the same 
resolution in fc-space by sampling a spectral distribution 
using the Box-Mueller method. The velocity field is trans¬ 
formed back to real-space using a Fast Fourier Transforma¬ 
tion (FFT), normalized such that the average velocity is of 
the desired Mach number. The velocities from the grid are 
transferred to the particle distribution using the nearest grid 
point (NGP) sampling kernel. 

To assess the impact of random motions near the res¬ 
olution scale, we need to measure the velocity power spec¬ 
trum within the SPH kernel. However, the accurate estima¬ 
tion of the velocity power spectrum of a particle distribu¬ 
tion close to the Nyqvist frequency is non-trivial, because of 
aliasing of the velocity power by the binning kernel (see dis¬ 
cussions in Jing 2005; Jasche et al. 2009; Cui et al. 2008). 
This can be compared to a problem in signal processing, 
where SPH represents an analogue signal and a grid a dig¬ 
ital signal representation of it. Aliasing is strongest at the 
smallest scales/largest modes, where the velocity power on 
the particles is modified by the shape of the binning kernel in 
configuration-space. To understand this effect and compare 
binning kernels, we take initial conditions with a full Kol¬ 
mogorov power spectrum (Pk oc k ~ n / 3 ) without performing 
a simulation and directly bin the velocity back to a grid 
using different kernels. After a forward fast Fourier trans¬ 


formation we radially average the velocity power in fc-space 
in 32 logarithmic bins. 

Fig. 10 shows the resulting power-spectra where the 
black line represents the original power spectrum. We show 
the kernels: Nearest Grid Point (NGP, red line), Cloud in 
Cell (CIC, green line), Triangular Shaped Cloud (TSC, vi¬ 
olet line), Daubechies scaling function of 20th order (D20, 
orange line), and the WC4 SPH kernel with 200 neighbors 
(SPH, brown line) (Hockney & Eastwood 1988; Daubechies 
1992; Dehnen & Aly 2012). We also show the NGP with 
two times oversampling (blue line), which was used in 
Bauer & Springel (2012). As vertical lines we show the 
Nyqvist frequency fcNyq = fVfc m i n , the WC4 smoothing scale 
k a and the WC4 kernel compact support fch sm i = 7t//i bm i 
(Dehnen & Aly 2012). 

During the binning process, the SPH kernel function 
conserves density to machine precision but not energy, i.e. 
binning with the SPH kernel is a diffusive process. The other 
kernels behave opposite, they conserve mass, scalar velocity 
and energy to less than one per cent but not SPH density 
and volume. Fig. 10 clearly shows that the D20 wavelet ker¬ 
nel minimizes aliasing for sufficiently homogeneous particle 
distributions (Cui et al. 2008). Our comparison also resolves 
the differences found in Bauer & Springel (2012) and Price 
(2012b), who use the twice oversampled NGP kernel and 
the standard SPH kernel, respectively. Prior studies based 
on the NGP kernel binning over-estimated the SPH noise, 
while SPH kernel based binning suppressed the real noise by 
aliasing. We conclude that all kernels except the D20 show 
substantial aliasing and it seems hard to draw definitive con¬ 
clusions from simulation results under this condition. Thus 
we make it our fiducial choice for this study. We note that in 
the presence of strong gradients in density the SPH kernel 
remains the only viable choice to obtain binned quantities, 
because it is the only kernel in our comparison that guar¬ 
antees a non-negative non-zero density in the entire simula¬ 
tion at all grid resolutions. This works reasonably well for 
a physical interpretation of velocity power-spectra, because 
motions below the SPH smoothing scale are caused by nu¬ 
merical effects (Price 2012b). 

3.6.2 Spectral evolution of turbulence 

To compare the ‘standard’ and the ‘new’ schemes we con¬ 
sider decaying turbulence within a periodic box. We seed 
compressive and solenoidal modes in the range of k £ 
[1.6,3.1], to obtain initial conditions appropriate for the 
galaxy and cluster environment, where turbulence is injected 
by merger infall on the scale of the halo core radius. We 
normalize the velocity fluctuations in the box such that the 
average velocity equals a Mach number of M = 0.1 and we 
do not time-average spectra. 

Fig. 11 (left and middle panels) shows the time evolu¬ 
tion of velocity power-spectra for the ‘standard’ scheme (left 
panel) and the ‘new’ scheme (middle panel). Here we also 
show the scale of the SPH kernel compact support (black 
vertical line) and the kernel smoothing scale (dotted ver¬ 
tical line). In-line with previous studies (Bauer & Springel 
2012 , their fig. 12), the ‘standard’ scheme does not develop 
a turbulent cascade and damps kinetic energy very quickly. 
Our ‘new’ scheme develops a cascade at large scales (small 
k) but then shows a depression of kinetic energy close to 
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Figure 9. Kelvin-Helmholtz instability. We show thin projections of specific entropy through the centre of the computational domain 
at times t = 3.0 and 6.0 (dynamical timescale /. kh ~ 3.4). In the ‘standard’ scheme, numerical surface tension as well as too much AV 
prevent the instability to develop and lead to an unphysical laminar behaviour of the fluid. In the ‘new’ scheme, our formulations of AV 
and AC promote the formation of roll-ups and onset of instability, while at late stages diffusion is dominating. 



Figure 10. Comparison of different particle to grid binning meth¬ 
ods acting on the same particle distribution sampling a Kol¬ 
mogorov velocity power spectrum in three dimensions (black line). 
We show the NGP kernel (red line), the twice oversampled NGP 
(blue line), the CIC (green line), the TSC (violet line), the D20 
(orange line) and the SPH WC4 (brown line). The vertical lines in¬ 
dicate the wave numbers corresponding to the Nyqvist frequency 
(solid line), the WC4 compact support (dashed line) and the WC4 
smoothing scale (dotted line). 


the kernel scale. This, again, is in-line with prior studies 
(Hopkins 2013). The damping of the spectrum at the later 
times appears self-similar, i.e. the shape of the spectrum 
does not change as energy decreases. Inside the kernel the 
typical build-up of thermal motions around the smoothing 


scale k a can be observed, but scales outside the kernel are 
not affected. 

In order to understand if the cause of the velocity de¬ 
pression at k ~ 10 is caused by the formulation of AV we per¬ 
form a test-run without any viscosity (Fig. 11, right panel). 
Throughout the whole evolution, the spectrum at the small¬ 
est k follows the Kolmogorov scaling, as expected. Once 
the turbulent cascade is established, the spectrum turns 
over at increasingly smaller scales, which is equivalent to an 
isotropization of kinetic energy inside the kernel and sub¬ 
sequent filling of larger scales with isotropic motions. This 
follows from the fact that SPH particles are subject to the 
pair-wise repulsive force (Price 2012a), and hence behave 
like a thermal gas below the kernel scale. Eventually, such a 
system will show a flat power spectrum as expected from the 
second law of thermodynamics. In our simulation, we also 
observe the depression outside the smoothing scale, which 
seems to be an intrinsic feature of SPH related to an energy 
transfer from outside the kernel to smaller scales, and not 
related to our formulation of AV. 

In our ‘new’ scheme including AV, thermal motions are 
well controlled inside the kernel scale, commonly referred to 
as ’kernel noise’. We argue that these motions are not spu¬ 
rious, because no additional energy is retained in them, as 
SPH is fully conservative and the spectrum decays roughly 
in a self-similar manner. If we define the kernel smoothing 
scale ( k a ) as the dissipation scale intrinsic to SPH, the ‘new’ 
scheme does not show a bottle-neck effect as found in grid 
codes at adjacent larger scales but a depression, roughly in 
the same range in k. This is in line with the results shown by 
Hopkins (2015), whose code uses a Riemann solver to formu¬ 
late noise-free AV on scales of the inter-particle separation 
to obtain grid-code behaviour. We note that the difference 
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Figure 11. Decaying subsonic turbulence. We show the build-up and decay of velocity power spectra for different schemes. The colours 
illustrate the time evolution of the spectra (sound-crossing time of about t 3 = 7.0. We initially distribute energy on the largest modes, 
which then develops a spectral distribution. In the ‘standard’ scheme (left panel) turbulent motions are almost completely suppressed by 
destructive impact of AV. In the ‘new’ scheme (middle panel) turbulent motions develop and a turbulent cascade is built. The diffusive 
character of AV is significantly changed and the velocity field as well as the kinetic energy are preserved. The spectra are then compared 
to a simulation without viscosity (right panel). 


in dissipation scale (d m in versus (Tkemei) translates into more 
resolution elements required by SPH compared to Eulerian 
methods, i.e. slower convergence. 

Fig. 12 shows thin slices through the centre of the sim¬ 
ulation box after one sound-crossing time. We visualize gas 
density (left panel) and vorticity (right panel). It can be 
clearly seen that our ‘new’ scheme resolves compressive and 
shearing velocity motions better than the ‘standard’ scheme 
The high-order derivatives of velocity lead to a more ac¬ 
curate estimation of vorticity and thus, limit the impact 
of AV and preserve kinetic energy and turbulent motions. 
Furthermore, the difference in AV and velocity dissipation 
between the two schemes becomes strikingly evident in the 
time evolution of kinetic energy (Fig. 13). The ‘standard’ 
scheme dissipates 90% of the kinetic energy budget within 
four sound-crossing times, while the ‘new’ scheme preserves 
energy better by a factor of 5. 

We conclude that our code performs comparably to 
modern implementations of SPH (Hopkins 2013; Price 
2012 a), even in the case of non-isothermal compressive 
and solenoidal decaying turbulence found in cosmologi¬ 
cal simulations. We show that the disagreement between 
Bauer & Springel (2012) and Price (2012a) is largely caused 
by technical differences, to solve it we propose a solution 
based on the D20 binning kernel. We also show that the 
downturn in the velocity power spectrum is not caused by 
the AV implementation. 


4 HYDRODYNAMICAL TESTS WITH 
GRAVITY 

We continue to evaluate the performance and accuracy of the 
two different SPH implementations with a second set of stan¬ 
dard problems. These second tests include hydrodynamical 
as well as gravitational forces and also take a cosmological 
time integration into account. 


4.1 Hydrostatic sphere 

We consider a sphere in hydrostatic equilibrium to study the 
SPH behaviour in combination with gravity in an ideally sta¬ 
ble system. We set up 88088 dark matter particles with indi¬ 
vidual masses of 2 ■ 10 9 Mq and 95156 gas particles with indi¬ 
vidual masses of 4.75-10 8 Mq. The total mass of the sphere is 

2.2 • 10 14 Mq and we use vacuum boundary conditions and a 
gravitational softening length of 12 kpc. We set up the initial 
equilibrium conditions following Komatsu & Seljak (2001), 
as described in Viola et al. (2008). We evolve the sphere adi- 
abatically and do not include cooling or heating mechanisms 
in this test. 

Fig. 14 shows the results of the test problem at vari¬ 
ous times. At first, the initial set-up (dotted lines) of the 
sphere is not yet in hydrostatic equilibrium and requires 
some time to settle. Once hydrostatic equilibrium is reached 
around time t = 2.6 (dashed lines for the ‘new’ scheme) 
all hydrodynamical schemes must preserve the structure of 
the sphere and the radial profiles towards the final sim¬ 
ulation time t = 7.6. Between all schemes, the thermal 
pressure profiles are indistinguishable balancing the gravi¬ 
tational pressure. However, the composition of the thermal 
pressure P = (7 — l)pit changes and the radial profiles for 
density p and internal energy u change significantly. The 
‘standard’ scheme (blue lines) reaches the highest central 
density and also features lowest central internal energy and 
entropy. However, the internal energy drops towards the cen¬ 
tre and no stable solution is reached at all. We suggest this 
behaviour to be a joint impact of pairing instability caused 
by the cubic spline kernel function and lack of internal en¬ 
ergy mixing. When introducing only AC, no gravitational 
limiter and no further SPH developments (green lines), the 
results marginally improve. Internal energy still drops to¬ 
wards the centre but this time because too much AC was 
introduced. In principle, AC leads to entropy cores but with¬ 
out the gravitational limiter over-mixing occurs and internal 
energy is transported from the centre to the outskirts along 
the pressure gradient. The addition of the gravitational lim¬ 
iter (brown lines) improves the results significantly. The di- 
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Figure 12. Decaying subsonic turbulence. We show thin slices through the centre of the simulation box after one sound-crossing time 
for both schemes. The panels correspond to gas density (left panel) and vorticity (right panel). The velocity field shows well developed 
turbulence consisting of compressive and shearing motions. The ‘new’ scheme is able to more accurately compute vorticity and suppress 
AV with the high-order Balsara limiter. 


vergence of profiles towards the centre is removed and a 
stable core of internal energy and entropy is reached. Fur¬ 
thermore, no numerically induced transport of heat takes 
place. Additionally, we perform a test run with introduc¬ 
ing only the WC4 kernel and no further SPH improvements 
(purple lines). In this run, the central density is lower than 
in the ‘standard’ scheme, most probably because the clump¬ 


ing of particles in the centre and the occurance of the pairing 
instability is suppressed by the WC4 kernel in contrast to 
the run with the cubic spline. This also leads to a plateau 
in internal energy in the centre of the sphere. At last, we 
show the results of the entire ‘new’ scheme (red lines) and 
find the results to remain stable with all additional modifi¬ 
cations to WC4 kernel, AV and time-step limiter. We find 
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Figure 13. Decaying subsonic turbulence. We show the total 
kinetic energy in the simulation box over time. 


this run to give the most stable radial profiles in time. This 
test confirms the importance of kernel functions immune to 
pairing instability and a proper implementation of AC and 
clearly shows the effects of particle clumping and over- and 
under-mixing in gravitationally virialized systems, which are 
important in cosmological simulations. 


an initial linear growth phase, the one-dimensional pertur¬ 
bation collapses and several strong shocks develop. Conve¬ 
niently, the Zel’dovich pancake has an analytical solution 
describing the evolution well up to the collapse, which we 
used to create the initial conditions of the simulation. The 
comoving position x of an initially unperturbated coordinate 
q at redshift z is given by 


x(q,z) = q- 


1 + z c sir i(kq) 
1 + z k 


(33) 


where k = 2-7r/A is the wavenumber of the perturbation with 
a wavelength of A. We numerically invert Eq. (33) to obtain 
q(x). The peculiar velocity corresponding to the initial dis¬ 
placement is given by 


c (a :,z) = -H 0 - 


1 + Z c 


sin (kq) 


(1 + *) 1 / 2 k ’ 

and the comoving density is given by 


p(x,z) = 


Po 


1 - TET cos(kq) ’ 


(34) 


(35) 


where po is the critical density, Ho the present-day Hubble 
constant and z c the redshift of collapse. Furthermore, the 
temperature evolves adiabatically up to the collapse as 


T(x,z ) = Ti 


Vi+*' 

\ 3 p(x,z) 

V 1 + Z i 

) po 


(36) 


4.2 Evrard collapse 

We consider the Evrard collapse (Evrard 1988) to study the 
SPH behaviour in the presence of dynamically important 
gravitational forces and collapse of gas. We initialize a sphere 
of gas with mass M = 1, radius R = 1 and density profile 
of p ~ r for r < R and use vacuum boundary conditions 
and a gravitational softening length of 0.005. We do not use 
an external gravitational potential, dark matter particles or 
radiative cooling and thus the cloud only self-gravitates on 
the free-fall time-scale. The gas is initially at rest and the 
thermal energy budget is orders of magnitude smaller than 
the gravitational binding energy. 

Fig. 15 shows the results of the test problem at time 
t = 0.8 We compare the SPH results to a reference solu¬ 
tion similar to Steinmetz & Mueller (1993). In density (left 
panel) as well as in velocity (middle panel) all schemes show 
similar trends. The general structure of the test is well re¬ 
produced; however, the shock front is slightly smoothed and 
broadened. The solution in density deviates slightly from 
the reference solution at the centre. The most striking dif¬ 
ferences can be seen in pre-shock entropy (right panel). The 
‘new’ scheme (red line) produces higher levels of entropy 
compared to the ‘standard’ scheme (blue line) We investi¬ 
gate this trend and perform another test simulation (purple 
line) with the ‘new’ scheme but the cubic spline kernel. This 
run is closest to the pre-shock reference solution. 

4.3 Zel’dovich pancake 

We consider the Zel’dovich pancake (Zel’dovich 1970) to 
study the SPH behaviour for the cosmological time inte¬ 
gration with Hubble function H(t) instead of time t. This 
test describes the evolution of a sinusoidal cosmological per¬ 
turbation in an expanding Einstein-de-Sitter universe. After 


where Zi is the initial redshift. We follow Bryan et al. (1995), 
Trac & Pen (2004) and Springel (2010b) in our test set-up 
and choose A = 64 Mpc/h, z c = 1, Zi = 100 and Ti = 100 K. 
In a fully three-dimensional box, we set up 256 3 dark matter 
particles of equal masses as well as 256 3 gas particles of equal 
masses. 

Fig. 16 shows the results of the test problem at two 
redshifts. In the top row we show the evolution of the pan¬ 
cake before the collapse at redshift a = 3.6, while it is still 
in the linear phase. The ‘standard’ (blue lines) and ‘new’ 
(red lines) schemes give comparable results in density con¬ 
trast (left column), temperature (middle column) and veloc¬ 
ity (right column). The simulated evolution agrees well with 
the analytical solution (green lines), which describes the lin¬ 
earized evolution of initial sinusoidal perturbation. At this 
stage of the collapse the test is dominated by gravitational 
forces and hydrodynamical forces are negligible. Therefore, 
we do not expect striking differences to arise between both 
schemes. Both capture the linear growth and adiabatic evo¬ 
lution well. 

In the bottom row we show the pancake at the final 
redshift 2 = 0. Again, we compare to the analytical solu¬ 
tion in the regions outside the central shock. In general, 
both schemes agree but we note the following differences. 
The peak density contrast is marginally lower in the ‘new’ 
scheme because of additional smoothing introduced by AC. 
However, the evolution of density contrast in the low-density 
regions is described better by the ‘new’ scheme and is re¬ 
solved with less noise. We recall the Sod shock tube (Sec¬ 
tion 3.1) and the Sedov blast wave (Section 3.2) problems, 
where we also find more accurately resolved density fields 
and lower peak densities. Concerning temperature, we find 
the central shock to be slightly broader in the ‘new’ scheme, 
which is caused by two effects. Firstly, the higher amount 
of viscosity within the shock leads to an earlier heating of 
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Figure 14. Sphere in hydrostatic equilibrium. We show radial profiles of density (left panel), internal energy (middle panel) and entropy 
(right panel) at three different times. At first, the initial conditions (dotted lines) must settle into hydrostatic equilibrium (dashed lines), 
which then remains stable for an extended period of time (solid lines). The stability of the sphere is determined by the differences 
occurring between the settled (only shown for the ‘new’ scheme) and the evolved state (shown for all schemes). In the ‘standard’ scheme 
(blue lines), pairing instability caused by the cubic spline kernel and lack of mixing lead to an incorrect central solution in density (too 
high) and internal energy (too low). The addition of the WC4 kernel (purple lines) prevents the formation of particle clumps at the 
centre and AC promotes fluid mixing but full AC without gravity limiter (green lines) leads to a numerical transport of internal energy 
outwards. The gravity limiter treats this behaviour (brown lines). The ‘new’ scheme (red lines) with WC4 and AC and gravity limiter 
significantly improves the radial profiles in all physical quantities. 




Figure 15. Evrard collapse. We show radial profiles of density, velocity and entropy at time t = 0.8 and compare to a reference piecewise 
parabolic grid computation (green lines). In principle, all schemes (red, blue and purple lines) show similar characteristics, but they differ 
as follows. The pre-shock entropy level is significantly higher in the ‘new’ scheme (red lines), which we attribute to the WC4 kernel with 
a larger smoothing size. A comparison run with the all improvements, but no WC4 kernel (purple lines), shows a lower pre-shock entropy 
level. 


particles and thus broadens the shocks. Secondly, the time- 
step limiting particle wake-up scheme captures highly active 
particles before they penetrate into inactive regions, which 
leads to a better fluid sampling and also shock broaden¬ 
ing. This can also be noticed in velocity, where the profile 
is slightly smoothed in the central region. In summary, our 
‘new’ scheme gives reasonable results in this cosmological 
test problem and the differences between both schemes are 
very small at redshift z = 0. Therefore, our new implemen¬ 
tation is ready to be applied to Idealized astrophysical prob¬ 
lems. 


5 ASTROPHYSICAL APPLICATIONS 

We complete the evaluation of the performance and accu¬ 
racy of the two different SPH implementations in idealized 
simulations of galaxy and galaxy cluster formation. 

5.1 Idealized galaxy formation 

In order to check if the ‘new’ scheme is numerically stable 
when coupled with a simple effective description of the inter¬ 
stellar medium, we consider the formation of an isolated disk 
from a cooling gas cloud embedded within a rotating dark 
matter halo. This idealized application also includes a the 
prescription for cooling, supernova feedback and star forma¬ 
tion of (Springel & Hernquist 2003). We focus on the differ- 
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Figure 16. Zel’dovich pancake. We show the evolution of density contrast (left column), temperature (middle column) and velocity (right 
column). We show the state of the pancake at an intermediate redshift z = 3.6, while it is still in the linear regime before the collapse and 
at the final redshift z = 0, when it is evolved well into the non-linear regime. Both SPH schemes agree well with the analytical solution 
during the linear growth phase. At the final redshift, the ‘new’ scheme resolves well the density contrast and yields a broader shock and 
a slightly smoother velocity. 


ences between both hydrodynamical schemes and therefore, 
we do not consider a cosmological environment or more ad¬ 
vanced physical processes such as black holes, stellar evolu¬ 
tion or metals. Numerical comparison simulations (see e.g. 
Scannapieco et al. 2012) are a common tool to study the 
impact of numerical schemes and physical modules. 

Within a computational domain of roughly 1 Mpc 3 we 
set up 4041345 particles resembling a Milky Way-like dark 
matter halo with a total mass of 1.8 • 10 12 Mq. We include 
4466429 gas particles with a total mass of 2.2 • 10 11 Mq, 
which corresponds to a baryon fraction of approximately 
11 per cent. The gravitational softening length is 450 pc. 
Initially, the distribution of dark matter follows a Navarro- 
Frenk-White profile (Navarro et al. 1997) and subsequently, 
we add the gaseous component similar to the set-up of the 
hydrostatic test (Section 4.1). The only change is that here, 
we give the gas and dark matter a rotational velocity which 
peaks at 180 km s _1 . Obviously, the initial hydrostatic equi¬ 
librium is broken by the onset of the gas cooling and we 
follow the evolution of the cloud for 10 Gyrs. 

Fig. 17 visualizes the spatial distribution of stars, where 
the colours represent the age of stars (top panels), and the 
spatial distribution of star forming gas, where the colours 
represent the star formation rate (bottom panels) at time 


t = 7.5 Gyr. We use the ray-tracing program SPLOTCH 
(Dolag et al. 2008; Jin et al. 2010) to create the images and 
choose a linear colour bar for stellar age, where the red end 
of the colour bar corresponds to stars older than 3 Gyrs and 
the blue end to recently formed stars. We visualize the star 
formation rate since it traces the (cold) gas within the disk 
and use a linear colour bar, which ranges from the minimum 
to the maximum value of star formation rate. We choose 
identical plot settings for the ‘standard’ scheme (left panels) 
and the ‘new’ scheme (right panels), which show striking 
morphological differences as follows. 

In the ‘standard’ scheme, the galaxy shows a prominent 
bulge containing a large fraction of the stellar population. 
The entire galaxy appears more spheroidal with a dominant 
bulge and the stellar disk is not well pronounced. We find 
similar features in the distribution of star formation. The 
gas disk is asymmetric and only shows little spiral struc¬ 
ture in the face-on view. In the edge-on projection the disk 
shows a rolling pin morphology. Both disks are dominated 
by bulges, but in the ‘new’ scheme the bulge is significantly 
less dominant. The bulge contains a smaller fraction of the 
stellar population and might eventually be associated with 
an elliptical bar structure. More as well as younger stars 
are present within the disk. The gas disk is symmetric and 
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shows a defined spiral structure. At this stage of code test¬ 
ing it is difficult to track down the impact of individual code 
changes. However, we assume the most significant differences 
are caused as follows. In the ‘standard’ scheme large amounts 
of AV might lead to dissipation of kinetic energy, loss of ro¬ 
tational support and numerical angular momentum trans¬ 
port. Additionally, the mixing problem and its associated 
numerical surface tension tend to confine cold and dense gas 
blobs. In the ‘new’ scheme, significantly smaller amounts of 
AV are applied and rotational support can be provided. Fur¬ 
thermore, the inclusion of AC promotes gas mixing between 
hot and cold phases. 

We continue with a more quantitative comparison of 
both schemes, which confirms our previous findings. Fig. 18 
shows density in vertical and radial direction (p z and p T ) of 
the stellar component as well as the associated vertical ve¬ 
locity dispersion (cr z ) at times 2.5 (dashed lines), 5.0 (dotted 
lines) and 7.5 Gyrs (solid lines). As seen in p z (left panel), 
the ‘new’ scheme produces a thinner stellar disk and as seen 
in p r (middle panel) the disk also extends to significantly 
larger radii. This trend is confirmed by <r z (right panel), 
which for the ‘standard’ scheme truncates at smaller radii 
than for the ‘new’ scheme. 

Fig. 19 shows vertical and radial profiles of gas density 
(p z and p r ) as well as the vertical gas velocity dispersion (a z ) 
of the cold gas at times 2.5 (dashed lines), 5.0 (dotted lines) 
and 7.5 Gyrs (solid lines). We employ a temperature crite¬ 
rion of T < 10 5 K to distinguish between cold and hot gas. p r 
decreases towards the centre since the gas within the bulge 
is hot and exceeds our temperature threshold. In the ‘new’ 
scheme, the distribution of cold gas is slightly more extended 
in vertical as well as radial direction. However, a z indicates 
a colder gas disk. Most probably, these features are a result 
of less numerically induced AV, angular momentum trans¬ 
port and depression of rotational support. Furthermore, the 
inclusion of AC allows mixing between gas phases and pro¬ 
motes dissociation of cold structures. 

5.2 Santa Barbara Cluster 

We carried out the Santa Barbara galaxy cluster 
(Frenk et al. 1999), which is a common reference simula¬ 
tion for cosmological hydrodynamical simulation codes. Al¬ 
though no analytic solutions exists, the cluster has been sim¬ 
ulated with a large variety of different codes. The simulation 
describes the formation of a massive dark matter halo, with 
virial mass of 1.2 • 10 15 Mq and virial radius of 2.8 Mpc. It 
is evolved in an Einstein-de Sitter cold dark matter cosmol¬ 
ogy with parameters of Qm = 1.0, Ha = 0.0, and Ho = 50 
km s' 1 Mpc -1 . We choose an initial redshift of 2 = 20 with 
a perturbed distribution of 256 3 dark matter particles and 
256 3 gas particles, each of equal masses (see also Frenk et al. 
1999, for a detailed description of the initial conditions), and 
follow the formation until redshift z = 0. 

Fig. 20 shows thin slices of gas density (left panels), tem¬ 
perature (middle panels) and entropy (right panels) defining 
the thermodynamical state of the hot intracluster medium 
(ICM) at redshift z = 0. For the definition of entropy we 
use S = T/n e 2 ^ 3 , which is commonly used in X-ray studies 
of the ICM (e.g. Kravtsov & Borgani 2012). From the maps 
we note the following interesting features. 

The gas density (left panels) tends to be smoother in the 


‘new’ scheme, which is mostly due to the effect of AC, which 
introduces entropy mixing among neighbouring gas parti¬ 
cles. In contrast, the ‘standard’ scheme produces a clumpy 
distribution of gas with gas inhomogeneities associated to 
stripping from merging haloes and cold blobs. These struc¬ 
tures are persistent in the hot ICM mainly due to the lack of 
mixing. In turn, these ’features’ of the ‘standard’ scheme pre¬ 
vent an efficient action of hydrodynamical instabilities such 
as Rayleigh-Taylor and Kelvin-Helmholtz instabilities that 
are spuriously inhibited. Quite remarkably, the clumps are 
much less evident in the ‘new’ scheme, which also produces 
a lower value for the central gas density. In the temperatures 
slices (middle panels) it becomes clear that gas clumps in the 
‘standard’ scheme correspond to objects of low temperature. 
As expected, the effect of introducing AC is the reduction 
of the degree of ICM thermal complexity. 

However, in the ‘new’ scheme the bow-shock, which is 
induced by the infall of a large substructure onto the main 
halo is better defined than in the ‘standard’ scheme. The 
bow-shock is located on the right hand side of the main halo 
centre (see Fig. 20). In fact, the WC4 kernel of the ‘new’ 
scheme captures the entropy discontinuity associated to the 
shock better. Most importantly, we note from the entropy 
maps (right panels) that the entropy level in the innermost 
region of the clusters increases in the ‘new’ scheme. 

Fig. 21 shows radial profiles of gas density, temperature 
and entropy in the same units as shown in Fig. 20 at redshift 
z — 0 for both SPH schemes. Furthermore, for a compari¬ 
son, we also include results obtained with the MASCLET 
grid code (see Quilis 2004). Both schemes produce quite 
consistent results at relatively large radii (> 200kpc/h) but 
show striking differences in the innermost regions. The ‘new’ 
scheme predicts a flatter central gas density profile, also with 
no evidence for the inversion of the temperature gradient 
produced by the ‘standard’ scheme. Density and tempera¬ 
ture profiles for the ‘new’ scheme combine to produce a flat 
entropy core, which extends out to ~ lOOkpc/h. In the ‘stan¬ 
dard’ scheme, the persistence of cold and dense clumps in 
the cluster atmosphere causes their low-entropy gas to sink 
to the centre of the cluster, thereby causing the continuous 
decrease of entropy. In the ‘new’ scheme, AC promotes mix¬ 
ing of gas phases and helps to dissolve low-entropy gas blobs 
within the hot ICM atmosphere and causes a higher entropy 
level to be established. 

The results for the ‘new’ scheme are remarkably simi¬ 
lar to those of the MASCLET code and, more in general, 
reported by Frenk et al. (1999); Vazza (2011); Power et al. 
(2014) for Eulerian codes. We point out that such a close 
agreement has been obtained without any tuning aimed at 
producing the entropy core predicted by Eulerian codes in 
cosmological simulations for the formation of galaxy clus¬ 
ters. The choice of parameters for the ‘new’ SPH scheme 
was only aimed at preventing the limitations of ‘standard’ 
SPH in terms of the description of discontinuities and ef¬ 
ficiency to capture hydrodynamical instabilities. Note also 
that these results for the Santa Barbara cluster are in qual¬ 
itative agreement with those obtained in the hydrostatic 
sphere (see Section 4.1). The behaviour of both schemes are 
in the same direction, even if they are less evident than in 
this Santa Barbara cluster test, due to the lack of the hierar¬ 
chical process of structure formation within a cosmological 
environment. 
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Figure 17. Idealized galaxy formation. We show the spatial distribution of stars, where the colours visualize the age of stars (top panels) 
and the spatial distribution of gas, where the colours visualize the star formation rate (bottom panels) at time t = 7.5 Gyr. In the 
‘standard’ scheme (left panels) the distribution of star formation is clumpy and the object appears bulgy. In the ‘new’ scheme (right 
panels) the distribution of star formation as well as the stellar component are more extended and pronounced in a disk-like structure. 
Furthermore, the size of the bulge is smaller and the distribution of young stars is more extended. We show a more quantitative comparison 
in Figs. 18 and 19. 


Additionally, we analyze the simulation also at red- 
shifts a > 0. In general, the profiles of the high-redshift 
haloes show the same behaviour as their low-redshift coun¬ 
terparts, provided that we choose quiet and virialized ob¬ 
jects. Objects, which host dynamically important shocks or 
undergo merger events show altered radial profiles because 
the timing of the mergers depends slightly on the simula¬ 
tion scheme. Therefore, a sensible comparison of the en¬ 
tire redshift-evolution and behaviour of both schemes during 


these violent phases of structure formation is not possible 
and requires controlled experiments. 


6 SUMMARY AND CONCLUSIONS 

In this paper we presented a novel implementation of the 
SPH scheme in the GADGET code, which provides im¬ 
proved accuracy for simulations of galaxies and large-scale 
cosmic structures. Since the first development of SPH great 
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Figure 18. Idealized galaxy formation. We show the vertical ( z ) density profile (left panel), the radial density profile (middle panel) 
and vertical velocity dispersion (right panel) of the stellar component. For the radial plots we use cylindrical bins. In the ‘new’ scheme 
(red lines), the galactic disk is more defined, extended and colder. 





Figure 19. Idealized galaxy formation. We show the vertical (z) density profile (left panel), the radial density profile (middle panel) and 
vertical velocity dispersion (right panel) of the cold gas component. For the radial plots we use cylindrical bins. We use a temperature 
criterion of T < 10 5 K to select the cold gas. In the ‘new’ scheme (red lines), the galactic disk is more defined, extended and colder. 


advancements have been made to improve the reliability and 
stability of this hydrodynamical scheme and, in particular, 
much effort has been spent in a proper treatment of disconti¬ 
nuities. We implemented and improved several of these mod¬ 
ifications of SPH into the developer version of GADGET-3, 
and tested them against a number of standard hydrodynam¬ 
ical problems, as well as first simple astrophysical applica¬ 
tions. The main modifications (see also Table 1) of this ‘new’ 
scheme, when compared to the ‘standard’ (see e.g. Price 
2012a) formulation of SPH, can be summarised as follows. 

• Artificial viscosity (AV) is introduced for a proper de¬ 
scription of shocks. It prevents particle interpenetration into 
unshocked regions and provides a regularisation of the par¬ 
ticle field, which supports a proper sampling of the fluid. 
First spatially constant low-order formulations of AV intro¬ 
duce viscosity not only at shocks, but also within unshocked 
regions and shearing flows, thereby leading to an overly 
viscous behavior and a too fast dissipation of kinetic en¬ 
ergy. Most commonly, the so called Balsara switch (Balsara 


1995) is used to reduce viscosity in shear flows, while fur¬ 
ther attempts were made to reduce AV where it is un¬ 
wanted (Morris & Monaghan 1997; Dolag et al. 2005). Re¬ 
cently, modern formulations of AV (Cullen & Dehnen 2010; 
Hu et al. 2014) improved greatly on a correct detection of 
shocks and use high-order gradient estimators to calculate 
divergence and curl of velocity from the full velocity gradient 
matrix instead of the classical SPH estimators. This allows 
shear flow limiters, such as the Balsara one, to work more 
accurately and suppress AV outside shocks and in shearing 
flows. In this way, kinetic energy is better preserved, thus 
helping simulating turbulent flows or hydrodynamical insta¬ 
bilities with higher accuracy. In our ‘new’ scheme, we com¬ 
pute the velocity gradients from the full velocity gradient 
matrix instead of low-order classic kernel derivatives. 

• Artificial conductivity (AC) is introduced to provide a 
proper fluid description at contact discontinuities. In fact, in 
the density-entropy formulation a spurious surface tension 
arises at discontinuities, which also suppresses the formation 
of instabilities and prevents mixing of different fluid phases. 
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Figure 20. Santa Barbara Cluster. In boxes with 1 Mpc/h side length, we show this slices of gas density (left panels), temperature 
(middle panels) and entropy (right panels) at redshift z = 0 for the ‘standard’ scheme (top row) and the ‘new’ scheme (bottom row). 
In the ‘new’ scheme significantly less dense and cold gas blobs are present, as AC promotes fluid mixing and blob dissociation. This 
promotes a smoother distribution of higher temperatures and entropies at the halo centre, which reduces the ICM complexity. 
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Figure 21. Santa Barbara Cluster. We show radial profiles of gas density (left panel), temperature (middle panel) and entropy (right 
panel). In each panel, we compare the results of the ‘standard’ (blue lines) and the ‘new’ (red lines) scheme to a reference solution (black 
lines) obtained with a piecewise parabolic grid computation with the MASCLET code (see Quilis 2004). The ‘standard’ scheme does not 
produce an entropy core but a diverging profile towards the halo centre. In the ‘new’ scheme an entropy core as well as stable temperature 
and density profiles are reached, which are all in good agreement to the grid code computation. 
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AV 

AC 

WakeUp 

Grav. 

Cosmo 

Radiat. 

Main result 

Sod shock tube 

/ 

/ 

- 

- 

- 

- 

Smooth density field / Sharp shock front / No pressure-blip 

Sedov blast wave 

/ 

/ 

/ 

- 

- 

- 

Sharp shock front / Central temperature profile 

Keplerian ring 

/ (imp.) 

- 

- 

- 

- 

- 

Stability of ring / No angular momentum transport 

Cold blob test 

/ 

/ 

- 

- 

- 

- 

No surface tension / Mixing of gas phases 

KH instability 

/ (imp.) 

/ 

- 

- 

- 

- 

Growth of perturbation / Mixing of gas phases 

Turbulent boxes 

/ (imp.) 

- 

- 

- 

- 

- 

Steady-state spectrum / Preservation of kinetic energy 

Hydrostatic test 

- 

/ (imp.) 

- 

/ 

- 

- 

Stability in radial profiles of density and entropy 

Evrard collapse 

/ (imp.) 

/ (imp.) 

/ 

/ 

- 

- 

Radial profiles in density and entropy 

Zel’dovich pancake 

/ 

/ (imp.) 

/ 

/ 

/ 

- 

Smooth density field / Sharp shock front 

Idealized galaxy 

/ (imp.) 

/ (imp.) 

- 

/ 

- 

/ 

Extended stellar and gas disks 

SB cluster 

/ (imp.) 

/ (imp.) 

/ 

/ 

/ 

- 

Formation of entropy core / Dissociation of cold blobs 


Table 2. Overview of our test problems. For each test we note the relative importance of a standard method (/) or an improved method 
(/(imp.)) of artificial viscosity (AV), artificial conductivity (AC), time-step limiter (WakeUp) and the physical processes involved beyond 
pure hydrodynamics such as gravity (Grav.), cosmological time integration (Cosmo), radiative cooling, star formation and supernova 
feedback (Radiat.). 


Lately, AC (see e.g. Price 2008) or pressure-entropy formula¬ 
tions (Hopkins 2013; Saitoh & Makino 2013) were proposed 
to overcome these issues. AC is applied at contact disconti¬ 
nuities and promotes the transport of heat between particles. 
However, in the presence of gravitationally induced pressure 
or temperature gradients, common AC schemes might lead 
to unwanted transport opposite to the gravitational force. 
Therefore, numerical limiters are necessary to be included 
in AC. In our ‘new’ scheme, we include locally adaptive AC 
to transport heat and treat contact discontinuities in SPH 
and we limit the amount of AC by correcting for gravita¬ 
tionally induced pressure gradients. While we demonstrated 
that our AC model is quite efficient at reducing such a sur¬ 
face tension, admittely a small residual effect is still present, 
and could potentially impact the long-term stability by over¬ 
diffusion. 

• As for the choice of the interpolating kernel, the com¬ 
monly employed cubic spline function has been shown to 
become easily unstable, which leads to spurious pairs of 
particles, incorrect gradient estimators and, in general, a 
poor fluid sampling. Therefore, a change of the kernel func¬ 
tion is highly recommended, where commonly the Wendland 
kernels (Dehnen & Aly 2012) are now used. In our ‘new’ 
scheme, we employ the Wendland C 4 kernel function with 
200 neighbours instead of a cubic spline with 64 neighbours. 
We calculate the density in a classic fashion from the mass 
distribution of particles and also compute the hydrodynam- 
ical forces with the density-entropy formulation. 

• At last, within supersonic shocks highly dynamical and 
computationally active particles can penetrate into regions 
containing computationally inactive particles causing dis¬ 
tortions in the fluid sampling and incorrect hydrodynami- 
cal solutions. In our ‘new’ scheme we use a particle wake- 
up time-step limiting scheme (see Saitoh & Makino 2009; 
Pakmor et al. 2012) as a solution, so as to shorten the time- 
steps whenever necessary and allow particles to become ac¬ 
tive earlier. 

To highlight the improvements associated to this ad¬ 
vanced SPH implementation, we investigate both the new 
and the original scheme in a variety of hydrodynamical stan¬ 
dard tests, with and without gravity. Furthermore, we study 


the behavior in the cosmological problem of the formation 
of galaxy cluster and enable simple prescriptions for radia¬ 
tive cooling, supernova feedback and star formation for a 
test simulation of an isolated rotating disk galaxy. Table 
2 presents an overview of our test problems and shows if 
SPH modules are important in a standard (/) or improved 
(/(imp.)) with respect to GADGET-SPH without our mod¬ 
ifications. Furthermore, we list the probed physical features 
of each test. 

The inclusion of AC in SPH changes the thermodynam¬ 
ical evolution of density, internal energy and pressure. Ad¬ 
ditionally, physical conduction (see e.g. Arth et al. 2014) is 
also sometimes employed in cosmological SPH simulations 
to promote (an)isotropic heat transport, which also helps to 
overcome the limitations of ‘standard’ GADGET-SPH. The 
joint effect of artificial conduction, introduced for purely nu¬ 
merical reasons, and of physical conduction awaits further 
investigations. 

In summary, the ‘new’ GADGET-SPH scheme pre¬ 
sented here performs better than the ‘standard’ one in ev¬ 
ery single of our test simulations. Therefore, it provides a 
much improved numerical description for weakly collision¬ 
less plasmas in cosmological simulations down to galactic 
scales. We base our future simulations of galaxies and large- 
scale cosmic structures on this updated formulation of SPH. 
We will also carry out detailed studies of galactic magnetic 
fields (see Beck et al. 2013, 2014) and the ICM with this ad¬ 
vanced method. In view of these applications, it is important 
to verify how this new SPH implementation performs when 
compared to other variants of SPH and to Eulerian codes. 
To this purpose, this SPH implementation participated to 
the nIFTy cosmology comparison project (Sembolini et al. 
2015), which compares the performances of different hydro- 
dynamical codes in cosmological re-simulations of galaxy 
clusters. In that comparison project, our code is shown to 
agree very well to both Eulerian codes and modern SPH im¬ 
plementations on the radial profiles of gas density, tempera¬ 
ture and entropy. Given the improvements in the description 
of hydrodynamics provided by the new SPH implementation 
presented here, we regard it as the core of an efficient code 
for modern simulations of cosmic structure formation. 
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